\documentclass{article}

\usepackage{geometry}
\usepackage{ctex}
\usepackage{amsmath}
\usepackage{lmodern}
\usepackage{listings}
\usepackage{hyperref}
%\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}    % 添加matlab代码宏
\usepackage{textcomp} % 必须加上，否则报错
\usepackage[dvipsnames]{xcolor}  % 更全的色系\usepackage{xcolor}
\usepackage{caption}
\usepackage{makecell}
\usepackage{graphicx}
\setmainfont{TeX Gyre Termes} %新增的


\author{20308063 \quad 黎逸雄}
\title{期中计算机作业\footnote{本次作业的所有源码放在个人gitee下，欢迎进行校验 \\ \url{https://gitee.com/lyxichigoichie/computer_hw_for_random_signal.git}}}
\geometry{a4paper, scale=0.8}
\lstset{
  language = Matlab,  %代码语言使用的是matlab
  %rulesepcolor=\color{red!20!green!20!blue!20},%代码块边框为淡青色
  rulesepcolor= \color{black},             % 代码块边框颜色
  keywordstyle=\color{blue!90}\bfseries, %代码关键字的颜色为蓝色，粗体
    numbers=left, % 显示行号
    numberstyle=\tiny,    % 行号字体
   commentstyle=\color{red!10!green!70},    % 设置代码注释的颜色
  showstringspaces=false,%不显示代码字符串中间的空格标记
  stringstyle=\ttfamily, % 代码字符串的特殊格式
  breaklines=true, %对过长的代码自动换行
  extendedchars=false,  %解决代码跨页时，章节标题，页眉等汉字不显示的问题
  escapebegin=\begin{CJK*}{GBK}{hei},escapeend=\end{CJK*},      % 代码中出现中文必须加上，否则报错
  texcl=true,
  frame = single}


\begin{document}
\maketitle
\section{反函数产生随机序列、概率密度函数}
\subsection{计算理想均值、方差}
\noindent 现有概率密度
\begin{equation*}
    f(x) = \begin{cases}
        \frac{x}{2}  &  0\leq x \leq 2\\
        0            & else
    \end{cases}
\end{equation*}
则有理想均值
\begin{equation*}
    \begin{split}
        E(X) &= \int_{-\infty}^{+\infty} xf(x) dx = \int_{0}^{2} \frac{x^2}{2} dx\\
             &= \frac{1}{6}x^3|_{0}^{2} = \frac{4}{3}
    \end{split}
\end{equation*}
二阶矩
\begin{equation*}
    \begin{split}
        E(X^2) &= \int_{-\infty}^{+\infty} x^2f(x) dx = \int_0^2 \frac{1}{2}x^3 dx\\
               &= \frac{1}{8}x^4|_{0}^{2} = 2
    \end{split}
\end{equation*}
理想方差
\begin{equation*}
    \begin{split}
        D(X) = E(X^2) - [E(X)]^2 = 2 - (\frac{4}{3})^2 = \frac{2}{9}
    \end{split}
\end{equation*}

\subsection{生成随机序列}
\noindent X的分布函数为
\begin{equation*}
    F_X(x) = \int_{-\infty}^{x} f(x) dx = \begin{cases}
        0 & x\leq 0\\
        \frac{1}{4}x^2 & 0 < x \leq 2\\
        1  & x>2
    \end{cases}
\end{equation*}
反函数为
\begin{equation*}
    F^{-1}(y) = \begin{cases}
        2\sqrt{y}  &  0 \leq y \leq 1\\
        \mbox{无意义}  & else
    \end{cases}
\end{equation*}
则可生成$X$的随机序列
\begin{equation*}
    X = F^{-1}(R)
\end{equation*}
其中$R$为在$(0,1)$均匀分布的随机变量

\begin{lstlisting}[caption = 产生均匀分布随机数并计算均值、方差]
clear all
clc
#生成在(0,1)均匀分布的1000个随机数序列
R = unifrnd(0,1,[1,1000]);
#求取X序列
X = 2.*sqrt(R);
#求取X序列的均值和方差
mean_value = mean(X)
variance = var(X)
#设置理想的均值和方差
std_mean_val = 4/3;
std_variance = 2/9;
#计算序列均值和方差与理想值之间的误差
error_mean = mean_value - std_mean_val
error_variance = variance - std_variance

#生成在(0,1)均匀分布的100000个随机数序列，并重新进行上述步骤
clear all
clc
R = unifrnd(0,1,[1,100000]);
X = 2.*sqrt(R);
mean_value = mean(X)
variance = var(X)
std_mean_val = 4/3;
std_variance = 2/9;
error_mean = mean_value - std_mean_val
error_variance = variance - std_varianceclear
\end{lstlisting}
计算结果为
\begin{table}[htbp]
    \centering
    \begin{tabular}{c*{4}{p{15ex}<\centering}}
        \hline \hline
        序列样本个数  &  均值  &  均值误差  &  方差  &  方差误差\\
        \hline
        $1\times 10^3$ & 1.3364  &  0.0031 & 0.2376 & 0.0154 \\
        %\hline
        $1\times 10^5$ & 1.3326 & $-7.44\times 10^{-4}$ & 0.2217 & $-5.21\times 10^{-4}$\\
        \hline
    \end{tabular}
\end{table}

\newpage
\subsection{绘制直方图及概率密度函数}
\begin{lstlisting}[caption = 绘制直方图及概率密度函数]
clear all
clc
R = unifrnd(0,1,[1,100000]);
X = 2.*sqrt(R);
#绘制直方图
h = histogram(X);
title('直方图')
print -dpng histogram.png

#求概率密度并绘制概率密度曲线
[f,xi]=ksdensity(X);
plot(xi,f,'LineWidth',3.0)
title('概率密度')
print -dpng probability_density.pngclear 
\end{lstlisting}
\begin{center}
    \begin{figure}[htbp]
        \begin{minipage}{0.45\linewidth}
            \centering
            \includegraphics[width=\linewidth]{figures/histogram.png}
            \caption{直方图}
        \end{minipage}
        \begin{minipage}{0.45\linewidth}
            \centering
            \includegraphics[width=\linewidth]{figures/probability_density.png}
            \caption{概率密度曲线}
        \end{minipage}
    \end{figure}
\end{center}

\newpage
\section{相关序列及其功率谱}
%\subsection{平稳正态序列的均值及方差的确定}
\noindent 由于正态随机序列的自相关函数满足
\begin{equation*}
    R_{X}(m) = \frac{1}{1-0.64} 0.8^{|m|}
\end{equation*}
即其自相关函数只与时间差有关，因此，只要再假设该随机序列的均值与时间无关，
则该序列就是一个平稳正态序列。因此，假设
\begin{equation*}
    m_X^2(n) = 0
\end{equation*}
%则可求得序列的方差
%\begin{equation*}
%    \begin{split}
%        R_X(0) &= \sigma_X^2 + m_X^2\\
%        \sigma_X^2 &= R_X(0) - m_X^2\\
%                   &= \frac{1}{1-0.64} - 0\\
%                   &= \frac{25}{9}
%    \end{split}
%\end{equation*}
%因此序列的分布为
%\begin{equation*}
%    X(n) \sim N(0,\frac{25}{9})
%\end{equation*}

若想得到具有某种相关关系的正态随机序列，则可通过生成N维正态
变量获得，该N维正态变量的协方差矩阵的元素为该相关函数减去均值的乘积。
由于已经假定均值为0，因此均值的乘积项为0。也即
\begin{equation*}
    \pmb{C} = \begin{bmatrix}
        R_X(0) & R_X(1) & \cdots & R_X(N-1)\\
        R_X(1) & R_X(0) & \cdots & R_X(N-2)\\
        \vdots & \vdots & \ddots & \vdots \\
        R_X(N-1) & R_X(N-2) & \cdots & R_X(0)
    \end{bmatrix}
    = \frac{\sigma^2}{1-a^2}\begin{bmatrix}
        1 & a & \cdots & a^{N-1}\\
        a & 1 & \cdots & a^{N-2}\\
        \vdots & \vdots & \ddots & \vdots \\
        a^{N-1} & a^{N-2} & \cdots & 1
    \end{bmatrix}
\end{equation*}
有如下变换
\begin{equation*}
    \pmb{X} = \pmb{AU}+\pmb{\mu}
\end{equation*}
其中$\pmb{X}$为具有相关关系$R_X(m)$的N维正态变量，$\pmb{\mu}$为均值矢量，
$\pmb{U}$为标准正态随机矢量，$\pmb{A}$可由以下分解得到
\begin{equation*}
    \pmb{C} = \pmb{AA^T}
\end{equation*}
则该正态随机序列可通过以下代码获得
\begin{lstlisting}[caption = 获得满足一定相关关系的正态随机序列]
#设定相关的常量
sigma = 1;
a = 0.8;
#设置产生的序列元素个数为500个
N = 500;
#将C视作对角阵D，下三角矩阵L和上三角矩阵L'的相加
D = eye(N);
L = zeros(N);
for i = 2:500
    for j = 1:i-1
        L(i,j) = a.^(i-j);
    end
end
C = sigma^2/(1-a^2).*(D + L' +L);
#将协方差矩阵C进行分解
A = (chol(C))';

#获得标准正态随机矢量
U = normrnd(0,1,500,1);
#获得相关正态随机序列并打印出来
X = A*U;
plot(X);
print -dpng 报告/figures/Correlated_normal_random_sequence.png
\end{lstlisting}
得到的序列的波形为
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.65\linewidth]{figures/Correlated_normal_random_sequence.png}
    \caption{相关正态随机序列波形}
\end{figure}

\begin{lstlisting}[caption = 绘制相关函数曲线]
#X是上段代码获得的序列
#实际相关函数
R = xcorr(X,'coeff');
plot(R,'b','LineWidth',1.0);
hold on;
#理论相关函数
m = -500 : 500;
RX = 1/(1-0.64).*(0.8).^(abs(m));
plot(RX,'r','LineWidth',1.0);
title('相关函数')
legend('实际曲线','理论曲线')
print -dpng 报告/figures/correlation_function.png
hold off
\end{lstlisting}

\newpage
\begin{lstlisting}[caption = 绘制功率谱密度曲线]
#X是上上段代码获得的序列
#功率谱密度
window = hann(length(X));
[Pxx,w] = periodogram(X,window);
plot(w,Pxx);
title('功率谱密度')
print -dpng 报告/figures/power_spectral_density.png
\end{lstlisting}

\begin{figure}[htbp]
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/correlation_function.png}
        \caption{相关函数曲线}
    \end{minipage}
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/power_spectral_density.png}
        \caption{功率谱密度}
    \end{minipage}
\end{figure}



\newpage
\section{ARMA模型1}
\subsection{模拟产生X(n)的500个观测点，画出波形并估计均值方差}
\begin{lstlisting}[caption = 模拟产生X(n)的500个观测点]
clear
clc
#系统的初始输入，为平稳白噪声的一个值
X = [2.*randn(1,1)];
#模拟AR过程产生500个观测点
for n = 2:500
    X = [X, -0.8*X(n-1)+2.*randn(1,1)];
end
#画出该样本函数的波形
n = 1:500;
plot(n,X)
title('500个观测点的波形')
print -dpng  报告/figures/500_observation_point_waveform.png
#求序列的均值和方差
x_mean = mean(X)
x_variance = var(X)
\end{lstlisting}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.7\linewidth]{figures/500_observation_point_waveform.png}
    \caption{500个观测点的样本函数波形}
\end{figure}
\noindent X(n)的均值和方差
\begin{equation*}
    \widehat{m_X} = 0.0438  \qquad \widehat{\sigma^2_X} = 9.7191
\end{equation*}

\subsection{画出理论的自相关函数和功率谱}
\noindent 一阶AR过程的自相关函数
\begin{equation*}
    R_X(m) = \frac{\sigma^2a^{|m|}}{1-a^2}
\end{equation*}
功率谱
\begin{equation*}
    G_X(\omega) = \frac{\sigma^2}{|1-ae^{-j\omega}|} = \frac{\sigma^2}{1+a^2-2a\cos\omega}
\end{equation*}
此处取$a=-0.8$

\begin{lstlisting}[caption = 画出理论自相关函数曲线]
a=-0.8;
m = -50:0.1:50;
R = 4.*a.^abs(m)/(1-a.^2);
plot(m,R)
title('理论自相关函数')
print -dpng 报告/figures/Theoretical_autocorrelation_function.png
\end{lstlisting}

\begin{lstlisting}[caption = 画出理论功率谱曲线]
a = -0.8;
omega = -6*pi:0.1:6*pi;
G = 4./(1+a.^2-2*a.*cos(omega));
plot(omega,G);
title('理论功率谱')
print -dpng 报告/figures/Theoretical_power_spectrum.png
\end{lstlisting}

\begin{figure}[htbp]
    \centering
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/Theoretical_autocorrelation_function.png}
        \caption{理论自相关函数曲线}
    \end{minipage}
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/Theoretical_power_spectrum.png}
        \caption{理论功率谱曲线}
    \end{minipage}
\end{figure}

\subsection{绘制实际自相关函数曲线及功率谱曲线}
\begin{lstlisting}[caption = 绘制实际自相关函数曲线]
#实际自相关函数
r = xcorr(X,'coeff');
plot(r);
title('实际自相关函数')
print -dpng 报告/figures/autocorrelation_function_for_three.png
\end{lstlisting}

\begin{lstlisting}[caption = 绘制实际功率谱密度曲线]
window = hann(length(X));
[Pxx,w] = periodogram(X,window);
plot(w,Pxx);
title('实际功率谱密度')
print -dpng 报告/figures/power_spectrum_for_three.png
\end{lstlisting}
\begin{figure}[htbp]
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/autocorrelation_function_for_three.png}
        \caption{实际自相关函数}
    \end{minipage}
    \begin{minipage}{0.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/power_spectrum_for_three.png}
        \caption{实际功率谱密度}
    \end{minipage}
\end{figure}


\end{document}